function [] = T_Q_sepdownload(startdate,enddate,uselocal,tqplot,useparfor)

%master script to load data and calculate T and Q
%Josh Crozier 2019

addpath('GISMO-master')
startup_GISMO

%may need to add javaclasspath.txt file to matlab pref folder 
%(find with prefdir, adjust to your directory)

%get data from sac files on local drive
if ~exist('uselocal','var')
    pstruct.uselocal = false; 
elseif isempty(uselocal)
    pstruct.uselocal = true;
else
    pstruct.uselocal = uselocal;
end
pstruct.localdrive = 'E:/';
% pstruct.localdrive = 'DawsonData/';

%plot tq and scaleogram
if ~exist('tqplot','var')
    pstruct.tqplot = true; 
elseif isempty(tqplot)
    pstruct.tqplot = false;
else
    pstruct.tqplot = tqplot;
end
%use parfor
if ~exist('useparfor','var')
    pstruct.useparfor = false; 
elseif isempty(useparfor)
    pstruct.useparfor = true;
else
    pstruct.useparfor = useparfor;
end

if pstruct.useparfor
    %setup parpool
    numcores = min([feature('numcores'),4]);
    %downloading is limiting factor, so no benifit to having more than ~3
    %workers
    poolobj = gcp('nocreate');
    if isempty(poolobj)
        newpool = true;
    elseif poolobj.NumWorkers ~= numcores
        newpool = true;
    else
        newpool = false;
    end
    if newpool
        poolobj = parpool(numcores);%,'AttachedFiles',{'IRIS-WS-2.0.18.jar'});
    end
end

if pstruct.useparfor
    pstruct.events_per_file = numcores-1; %number of events to save at a time
else
    pstruct.events_per_file = 1; %number of events to save at a time
end

main_starttime = datetime;


%% setup to download and process data
%% general settings
pstruct.cwt_save = false; %save cwt data
pstruct.fmplot = false; %plot first motions
pstruct.cwt_compress = {seconds(60),200}; % output compressed cwt res {timewindow,numfreqbins} (empty for no compression)
pstruct.gapdur = seconds(4); %allowable data gap (s)

pstruct.water_level = 1E1; %max magnification relative to min magnification
%reccommend around 10 for HVO broadband network

if ~exist('startdate','var')
    if ~pstruct.uselocal
        pstruct.startdate = datetime(2012,1,1,0,0,0);
    else
        pstruct.startdate = datetime(2008,1,1,0,0,0);
    end
    'using different startdate'
    pstruct.startdate = datetime(2018,5,6,15,30,0);
elseif isempty(startdate)
    if ~pstruct.uselocal
        pstruct.startdate = datetime(2012,1,1,0,0,0);
    else
        pstruct.startdate = datetime(2008,1,1,0,0,0);
    end
else
    pstruct.startdate = startdate;
end


if ~exist('enddate','var')
    if ~pstruct.uselocal
        pstruct.enddate = datetime(2019,1,1,0,0,0);
    else
        pstruct.enddate = datetime(2012,1,1,0,0,0);
    end
    'using different enddate'
   pstruct.enddate = pstruct.startdate + minutes(30);
elseif isempty(enddate)
    if ~pstruct.uselocal
        pstruct.enddate = datetime(2019,1,1,0,0,0);
    else
        pstruct.enddate = datetime(2011,12,29,0,0,0);
    end
else
    pstruct.enddate = enddate;
end


% '!!!test prefix!!!'
pstruct.foldername = strcat('TQnew_',datestr(pstruct.startdate,30),'_',datestr(pstruct.enddate,30)); %filename for output
pstruct.TQ_table_foldername = strcat(pstruct.foldername,'/TQ_table'); %filename for output
pstruct.cwt_foldername = strcat(pstruct.foldername,'/cwt'); %filename for output

pstruct.window = minutes(30*1); %length of data to examine [s] 
%(extra padding will actually need to be loaded also, done below)
%values of 12hr or larger will cause problems if reading local files

pstruct.networks = "HV";
if pstruct.uselocal
    pstruct.stationlist = ["NPB","SRM","NPT","OBL","WRM","SDH","UWE","UWB",...
        "SBL","KKO","RIMD"];
else
    pstruct.stationlist = ["NPT","OBL","WRM","SDH","UWE","UWB",...
        "SBL","KKO","RIMD"];
end

% 'using limited stations'
% pstruct.stationlist = ["NPT"];

pstruct.locations = "*";
pstruct.channellist = ["HHE","HHN","HHZ"];


%% using synthetics
pstruct.use_synthetics = false;
if pstruct.use_synthetics
    disp('!!!Using Synthetics!!!')
    pstruct.synthetic_T = seconds([20]); %array of oscillation periods
    pstruct.synthetic_Q = [20]; %array of oscillation quality-factors
    pstruct.synthetic_A = 1E3*[1]; %array of oscillation amplitudes
    pstruct.synthetic_onset = seconds(60*[5]); %array of oscillation onset times
    pstruct.synthetic_step = 1E3*[2.5];%step amplitude
    pstruct.synthetic_noise = 0; %white noise amplitude (ratio of total signal amp)
    pstruct.synthetic_sourcenoise = 0; %white noise amplitude in sourcetimefinc (ratio of total signal amp)
    pstruct.synthetic_mu = 10E9; %shear modulus
    pstruct.synthetic_nu = 0.25; %poisson
    pstruct.synthetic_lat = 19.405402; %latitude
    pstruct.synthetic_lon = -155.281041; %longitude
    pstruct.synthetic_elevation = 0; %elevation (asl)
    pstruct.cwt_save = false;
    pstruct.events_per_file = 1;
    pstruct.tqplot = true;
    pstruct.fmplot = true;
end


%% Problem dependent TQ calc settings
pstruct.Tmax = seconds(60); % max period
pstruct.Tmin = seconds(2); % min period
pstruct.minSTALTA = 3.0; %min amp/LTA ratio
pstruct.minampstd = 3.0; %min log amp of peaks (rel to std of log LTA)
pstruct.minpromstd = 0; %min prominence of peaks (rel to std of log LTA)
pstruct.minTseprat = 1.07; %minimum ratio of periods between peaks
pstruct.mintsep = seconds(200); %minimum separation between peaks in time
pstruct.decaythresh = 0.5; %min remaining amplitude fraction per cycle (controls min Q)
%will need to be increased with wider wavelets to avoid artificial picks
pstruct.minhalflife = seconds(6); %minumum oscillation half life 
pstruct.LTA_window = seconds(200); %windowspan for long-term avg

%these nest 2 parameters shouldn't be used with the current detection setup,
%which already chooses the peak with the highest amp integrated over 2
%cycles, but both params can be left at 1 just in case
pstruct.Aweight = 1; %weight of mode "goodness" on amplitude
pstruct.gweight = 1; %weight mode "goodness" on decay rate


%% Mostly general TQ calc settings (probably don't need much adjustment)
pstruct.VoicesPerOctave = 48; %for cwt, increase for better freq resolution
pstruct.fitstartoffsetcycles = 1.0;%offset to start expfit,recommend 0-2 cycles
pstruct.LTAcycleoffset = 4.0; %amount to shift LTA (accouns for transform delay)
%will depend to some extent on wavelet properties
pstruct.WaveletWeights = [1,0,0,1]; %[morse, bump, spec, narrowmorse], weights on different transforms
pstruct.QWaveletWeights = [0,0,0,1]; %[morse, bump, spec, narrowmorse], for expfit
% Morse (3-60) and amor similar but morse bit better
% bump yeilds better spectral res (even than morse 3-120),
% but bump has worse temporal res and some artificial banding
% spectrogram not as useful as either
pstruct.WaveletParameters = [3,120]; %morse wavelet params for cwt
% recommend [3, 60-120], 120 gives better f resolution but artificially
% increases Q for low Q events, so using additional narrow wavelet helps
% mitigate
pstruct.NarrowWaveletParameters = [3,60];
% 12-30 seems to be good balance, if too narrow will capture too many other
% frequencies to yeild robust estimates of Q
pstruct.phase_WaveletParameters = [3,60];
% for checking that phase is consistent across event, needs to be narrow
% enough to detect gaps but if too narrow will capture other frequencies
pstruct.specwindowtime = seconds(120); %spectrogram window (also governs edge trimming)
pstruct.fitcycles = 8; %sets number of cycles to exp fit (reccomend ~4-10)
pstruct.phasefitcycles = 6; %sets number of cycles to exp fit (reccomend ~4-10)
pstruct.expmodelvars = 'b_under'; %variables to fit in "(d-c)*exp(bx)+c"
% 'bcd' for full 3 param
% 'bc' for fixed amp variable offset
% just 'b' for fixed amp and offset
% 'bd' for fixed offset variable amp
% 'b_under' for fixed amp + offset and fit under all data points
% 'b_under' seems to be most robust option
pstruct.weight_decay_scale = 1; %(only used if not 'b_under') decaying exp fit weights, smaller for faster decay
% recommend around 1
pstruct.c_quantile = inf; %quantile at which to set constant offset for exp fits (for not 'b_under')
% recommend around 0.5-0.01
pstruct.c_std_frac = inf; %fraction of range beneath min to set c (for 'b_under')
% recommend inf, if inf then c = 0
pstruct.minadjrsquare = -Inf; %min adjusted r squared value for exp fit
% recommend -Inf for 'b_under' 
pstruct.rsquaredweight = 0; %weight mode "goodness" on exp fit r^2
% recommend 0 for 'b_under', small value relative to other weights or zero otherwise
pstruct.promweight = 0; %weight mode "goodness" on prominence
% recommend 0 for all cases
pstruct.ampcycles = 3; %number of cycles to sum over to check that is a local max
% recommend around 0.5-3
pstruct.fm_LTA_window_cycles = 3; %number cycles for STA/LTA first motion pick
pstruct.fm_minampfrac = 0.5; %min amp for first motion pick relative to max 
%(prevents picking artifacts)
pstruct.startbounds = pstruct.specwindowtime/2 + ...
    [pstruct.LTA_window, pstruct.Tmax*pstruct.fitcycles]; %peak start time bounds

pstruct.startshift = pstruct.startbounds(1);
pstruct.endshift = pstruct.window + pstruct.startbounds(2);
pstruct.taperfrac = seconds(pstruct.startshift/seconds(pstruct.window + pstruct.startshift))/2;

%% main
% stations = join(stationlist,",");
% channels = join(channellist,",");
% pstruct.useirisfetch = false; %false for GISMO instead of irisfetch
%currently only GISMO works

%table names for each channel
pstruct.channel_station_str = strings(1,length(pstruct.stationlist)*length(pstruct.channellist));
for statin = 1:length(pstruct.stationlist)
    for chanin = 1:length(pstruct.channellist)
        in = length(pstruct.channellist)*(statin - 1) + chanin;
        pstruct.channel_station_str(in) = strcat(pstruct.networks,'_',...
            pstruct.stationlist(statin),...
            '_', pstruct.channellist(chanin));
    end
end
pstruct.channel_station_list = pstruct.channel_station_str;

% if ~pstruct.useirisfetch
    %channel tags for GISMO
    pstruct.ctags = ChannelTag.empty(length(pstruct.stationlist)*length(pstruct.channellist),0);
    for statin = 1:length(pstruct.stationlist)
        for chanin = 1:length(pstruct.channellist)
            in = length(pstruct.channellist)*(statin - 1) + chanin;
            pstruct.ctags(in) = ChannelTag(pstruct.networks,pstruct.stationlist(statin), ...
                pstruct.locations, pstruct.channellist(chanin));
        end
    end
% end

% if useirisfetch
%     % make sure the java jar is in the path, this need only be done once per MATLAB session
%     %javaaddpath("IRIS-WS-2.0.18.jar");
% else
    if ~pstruct.uselocal
        %seismic data from iris
        pstruct.ds = datasource('irisdmcws');
        pstruct.localdrive = '';
        %local reponses (due to iris bug)
        pstruct.resps = {};
        for ii = 1:length(pstruct.ctags)
            temp_filename = strcat('Responses/',pstruct.ctags(ii).station,'-',...
                pstruct.ctags(ii).channel,'-sacpz_iris.txt');
            pstruct.resps{ii} = sacpz(fileread(temp_filename));
        end
    else
        pstruct.ds = '';
        %load from local drive
        %also load responses
        for ii = 1:length(pstruct.ctags)
            temp_filename = strcat('Responses/',pstruct.ctags(ii).station,'-',...
                pstruct.ctags(ii).channel,'-sacpz.txt');
            resp = sacpz(fileread(temp_filename));
            resp = struct('zeros',resp.z,...
                'poles',resp.p,...
                'constant',resp.k);
            if ii == 1
                pstruct.responses = resp;
            else
                pstruct.responses(ii) = resp;
            end
        end
    end
%     javaaddpath("IRIS-WS-2.0.18.jar");
%     antelopeds = '/aerun/sum/db/dbmaster/master_stations';

%create directories for saving output
if exist(pstruct.foldername,'file') ~= 7
    mkdir(pstruct.foldername)
end
if exist(pstruct.TQ_table_foldername,'file') ~= 7
    mkdir(pstruct.TQ_table_foldername)
end
if exist(pstruct.cwt_foldername,'file') ~= 7
    mkdir(pstruct.cwt_foldername)
end

%process batches of events at a time before saving
batchcount = ceil((pstruct.enddate - pstruct.startdate)/...
    (pstruct.window*pstruct.events_per_file));
cur_trace_cell = cell(pstruct.events_per_file+1,1);
next_trace_cellcell = cell(pstruct.events_per_file+1,1);
TT_datetime_local = [];
for batchin = 0:batchcount
    
    next_TT_datetime_local = (pstruct.events_per_file*pstruct.window*batchin + ...
        pstruct.window*(0:1:pstruct.events_per_file-1)) + pstruct.startdate;
    
    TQ_tab_cell = cell(pstruct.events_per_file+1,1);
    if pstruct.tqplot
        TQ_plot_cell = cell(pstruct.events_per_file+1,1);
    end
    
    cwt_cell = cell(pstruct.events_per_file + 1,1);
    
    % download and process
    if pstruct.useparfor
        parfor processin = 1:pstruct.events_per_file + 1
            if processin == pstruct.events_per_file + 1 && batchin < batchcount
                [next_trace_cell_temp] = ...
                    load_fun(pstruct,...
                        next_TT_datetime_local,...
                        main_starttime,...
                        batchin);
                next_trace_cellcell{processin} = next_trace_cell_temp;
                
            elseif batchin > 0
                [TQ_plot_temp, TQ_tab_temp,cwt_struct] = ...
                    process_fun(pstruct,...
                        cur_trace_cell{processin},...
                        TT_datetime_local,...
                        batchin,...
                        processin);
                TQ_plot_cell{processin} = TQ_plot_temp;
                TQ_tab_cell{processin} = TQ_tab_temp;
                cwt_cell{processin} = cwt_struct;
            end
        end %end loop over events
    
    else %if not useing parallell
        for processin = 1:pstruct.events_per_file + 1
            if processin == pstruct.events_per_file + 1 && batchin < batchcount
                [next_trace_cell_temp] = ...
                    load_fun(pstruct,...
                        next_TT_datetime_local,...
                        main_starttime,...
                        batchin);
                next_trace_cellcell{processin} = next_trace_cell_temp;
                
            elseif batchin > 0
                [TQ_plot_temp, TQ_tab_temp,cwt_struct] = ...
                    process_fun(pstruct,...
                        cur_trace_cell{processin},...
                        TT_datetime_local,...
                        batchin,...
                        processin);
                TQ_plot_cell{processin} = TQ_plot_temp;
                TQ_tab_cell{processin} = TQ_tab_temp;
                cwt_cell{processin} = cwt_struct;
            end
        end %end loop over events
    end %end if parfor
    
    if batchin < batchcount
        %update traces for next batch
        cur_trace_cell = cell(pstruct.events_per_file+1,1);
        for eventin = 1:pstruct.events_per_file
            cur_trace_cell{eventin} = next_trace_cellcell{pstruct.events_per_file+1}{eventin};
            next_trace_cellcell{pstruct.events_per_file+1}{eventin} = []; %saves memory
        end
        next_trace_cellcell = cell(pstruct.events_per_file+1,1);
        TT_datetime_local = next_TT_datetime_local;
    end

    if batchin > 0
        if pstruct.tqplot
            for jj = 1:pstruct.events_per_file
                if isgraphics(TQ_plot_cell{jj})
                    figure(TQ_plot_cell{jj})
                end
            end
        end

        %%
        %save
        validinds = false(size(TQ_tab_cell));
        for validin = 1:length(validinds)
            if ~isempty(TQ_tab_cell{validin})
                validinds(validin) = true;
            end
        end
        TQ_tab = cat(1, TQ_tab_cell{validinds});
        if ~isempty(TQ_tab)
            writetable(TQ_tab, strcat(pstruct.TQ_table_foldername,'/',num2str(batchin),'.csv'));
        end
        if pstruct.cwt_save&&~isempty(pstruct.cwt_compress)
            cwt_struct = [];
            for jj = 1:pstruct.events_per_file
                if isempty(cwt_struct)
                    if ~isempty(cwt_cell{jj})
                        cwt_struct = struct();
                        cwt_struct.time = cwt_cell{jj}.time;
                        cwt_struct.cwt = cwt_cell{jj}.wt;
                        cwt_struct.s = cwt_cell{jj}.s;
                        cwt_struct.numchannels = cwt_cell{jj}.numchannels...
                            *zeros(1,length(cwt_cell{jj}.time));
                    end
                else
                    if ~isempty(cwt_cell{jj})
                        cwt_struct.time = [cwt_struct.time, ...
                            cwt_cell{jj}.time];
                        cwt_struct.cwt = [cwt_struct.cwt,...
                            cwt_cell{jj}.wt];
                        cwt_struct.numchannels = [cwt_struct.numchannels,...
                            cwt_cell{jj}.numchannels+zeros(1,length(cwt_cell{jj}.time))];
                    end
                end
            end
            if ~isempty(cwt_struct)
                save(strcat(pstruct.cwt_foldername,'/cwt_',...
                    datestr(TT_datetime_local(1),30)), 'cwt_struct')
            end
        end
    end
         
end %end loop over files (event batches)


%% combine tables
% for batchin = 1:batchcount
%     if exist(strcat(pstruct.TQ_table_foldername,'/',num2str(batchin),'.csv'),'file') == 2
%         TQ_tab_local = readtable(strcat(pstruct.TQ_table_foldername,'/',num2str(batchin),'.csv'));
%         if batchin == 1
%             TQ_tab = TQ_tab_local;
%         else
%             TQ_tab = cat(1,TQ_tab,TQ_tab_local);
%         end
%     end
% end
% %save
% writetable(TQ_tab, strcat(pstruct.TQ_table_foldername,'/total','.csv'));

end %end main function


% %% for saving
% function cwt_savefun(filename, TQ_tab_local, wt, time, s, coi)
%     save(filename, 'TQ_tab_local', 'wt', 'time', 's', 'coi')
% end


